#include <life/options.hpp>
#include <life/lifecore/application.hpp>

#include <life/lifealg/backend.hpp>

#include <life/lifediscr/functionspace.hpp>

#include <life/lifepoly/im.hpp>

#include <life/lifefilters/gmsh.hpp>
#include <life/lifefilters/gmshtensorizeddomain.hpp>
#include <life/lifefilters/exporter.hpp>
#include <life/lifepoly/lagrange.hpp>
#include <life/lifepoly/crouzeixraviart.hpp>

#include <life/lifemesh/elements.hpp>

#include <life/lifevf/vf.hpp>

#include <life/lifefilters/geotool.hpp>

#include <life/lifediscr/operatorlinear.hpp>
#include <life/lifediscr/operatorinterpolation.hpp>

/**
 * This routine returns the list of options using the
 * boost::program_options library. The data returned is typically used
 * as an argument of a Life::Application subclass.
 *
 * \return the list of options
 */

namespace Life
{
namespace FBM_Cavity
{

inline
Life::po::options_description
makeOptions()
{
    Life::po::options_description stokesoptions("Stokes options");
    stokesoptions.add_options()
        ("penal", Life::po::value<double>()->default_value( 0.5 ), "penalisation parameter")
        ("f", Life::po::value<double>()->default_value( 0 ), "forcing term")
        ("mu", Life::po::value<double>()->default_value( 1.0 ), "reaction coefficient component")
        ("hsize", Life::po::value<double>()->default_value( 0.1 ), "first h value to start convergence")
        ("hlocsize", Life::po::value<double>()->default_value( 0.1 ), "first h value to start convergence")
        ("nbiter", Life::po::value<int>()->default_value( 50 ), "0 = strong Dirichlet, 1 = weak Dirichlet")
        ("bctype", Life::po::value<int>()->default_value( 0 ), "0 = strong Dirichlet, 1 = weak Dirichlet")
        ("bccoeff", Life::po::value<double>()->default_value( 100.0 ), "coeff for weak Dirichlet conditions")
        ("export-matlab", "export matrix and vectors in matlab" )
        ;
    return stokesoptions.add( Life::life_options() ) ;
}


/**
 * This routine defines some information about the application like
 * authors, version, or name of the application. The data returned is
 * typically used as an argument of a Life::Application subclass.
 *
 * \return some data about the application.
 */
inline
Life::AboutData
makeAbout()
{
    Life::AboutData about( "fbm_stokes_cavity" ,
                           "fbm_stokes_cavity" ,
                           "0.1",
                           "Stokes equation on simplices or simplex products",
                           Life::AboutData::License_GPL,
                           "Copyright (c) 2009 Universite de Grenoble 1 (Joseph Fourier)");

    about.addAuthor("Christophe Prud'homme", "developer", "christophe.prudhomme@ujf-grenoble.fr", "");
   return about;

}


/**
 * \class Stokes class
 * \brief solves the stokes equations
 *
 */
template<int Dim,
         typename BasisU,
         typename BasisP,
         template<uint16_type,uint16_type,uint16_type> class Entity>
class Stokes
    :
        public Application
{
    typedef Application super;
public:


    typedef double value_type;

    typedef Backend<value_type> backend_type;
    typedef boost::shared_ptr<backend_type> backend_ptrtype;

    /*matrix*/
    typedef typename backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
    typedef typename backend_type::vector_ptrtype vector_ptrtype;

    /*mesh*/
    typedef Entity<Dim,1,Dim> convex_type;
    typedef Mesh<convex_type> mesh_type;
    typedef boost::shared_ptr<mesh_type> mesh_ptrtype;

    /*basis*/
    typedef BasisU basis_u_type;
    typedef BasisP basis_p_type;
    typedef Lagrange<0, Scalar> basis_l_type;
    typedef bases<basis_u_type,basis_p_type, basis_l_type> basis_type;

    /*space*/
    typedef FunctionSpace<mesh_type, basis_type> space_type;
    typedef boost::shared_ptr<space_type> space_ptrtype;


    typedef typename space_type::template sub_functionspace<0>::type::element_type space_u_type;
    typedef boost::shared_ptr<space_u_type> space_u_ptrtype;
    typedef typename space_type::template sub_functionspace<1>::type::element_type space_p_type;
    typedef boost::shared_ptr<space_p_type> space_p_ptrtype;
    typedef typename space_type::template sub_functionspace<2>::type::element_type space_l_type;
    typedef boost::shared_ptr<space_l_type> space_l_ptrtype;

    //typedef FunctionSpace<mesh_type, bases<basis_p_type> > space_p_type;
    //typedef boost::shared_ptr<space_p_type> space_p_ptrtype;
    //typedef typename space_p_type::element_type element_p_type;

    BOOST_MPL_ASSERT( ( boost::is_same<typename space_type::bases_list, basis_type> ) );
    BOOST_MPL_ASSERT( ( boost::is_same<typename mpl::at<typename space_type::bases_list,mpl::int_<0> >::type, basis_u_type> ) );
    BOOST_MPL_ASSERT( ( boost::is_same<typename mpl::at<typename space_type::bases_list,mpl::int_<1> >::type, basis_p_type> ) );
    BOOST_MPL_ASSERT( ( boost::is_same<typename mpl::at<typename space_type::bases_list,mpl::int_<2> >::type, basis_l_type> ) );

    /* functions */
    typedef typename space_type::element_type element_type;
    typedef typename element_type::template sub_element<0>::type element_0_type;
    typedef typename element_type::template sub_element<1>::type element_1_type;
    typedef typename element_type::template sub_element<2>::type element_2_type;

    /* particule */
    typedef FunctionSpace<mesh_type, bases<Lagrange<0, Scalar> >, Discontinuous > p0_functionspace_type;
    typedef boost::shared_ptr<p0_functionspace_type> p0_functionspace_ptrtype;
    typedef typename p0_functionspace_type::element_type p0_element_type;

    /* interpolation */
    typedef OperatorInterpolation<  typename space_type::template sub_functionspace<0>::type::element_type ,
                                    typename space_type::template sub_functionspace<0>::type::element_type > opinterp_u_type;
    typedef OperatorInterpolation<  typename space_type::template sub_functionspace<1>::type::element_type ,
                                    typename space_type::template sub_functionspace<1>::type::element_type > opinterp_p_type;
    typedef OperatorInterpolation<  typename space_type::template sub_functionspace<2>::type::element_type ,
                                    typename space_type::template sub_functionspace<2>::type::element_type > opinterp_l_type;
    //typedef OperatorInterpolation<  typename space_type::template sub_functionspace<1>::type::element_type,//space_p_type,
    //                                space_p_type //p0_functionspace_type
    //> opinterp_p0_type;

    //typedef OperatorInterpolation<space_type,space_type> opinterp_type;
    typedef boost::shared_ptr<opinterp_u_type> opinterp_u_ptrtype;
    typedef boost::shared_ptr<opinterp_p_type> opinterp_p_ptrtype;
    typedef boost::shared_ptr<opinterp_l_type> opinterp_l_ptrtype;
    //typedef boost::shared_ptr<opinterp_p0_type> opinterp_p0_ptrtype;

    //typedef FsFunctionalLinear<typename space_type::template sub_functionspace<0>::type::element_type> flin_type;
    //typedef boost::shared_ptr<flin_type> flin_ptrtype;

    //typedef FsFunctionalLinear< space_p_type> flin_p0_type;
    //typedef FsFunctionalLinear<p0_functionspace_type> flin_p0_type;
    //typedef boost::shared_ptr<flin_p0_type> flin_p0_ptrtype;

    /* export */
    typedef Exporter<mesh_type> export_type;
    typedef boost::shared_ptr<export_type> export_ptrtype;

    Stokes( int argc, char** argv, AboutData const& ad, po::options_description const& od )
        :
        super( argc, argv, ad, od ),
        M_backend( backend_type::build( this->vm() ) ),
        M_backend_Loc( backend_type::build( this->vm() ) ),
        meshSize( this->vm()["hsize"].template as<double>() ),
        meshLocalSize( this->vm()["hlocsize"].template as<double>() ),
        exporter( Exporter<mesh_type>::New( this->vm(), this->about().appName() ) )
    {
        mu = this->vm()["mu"].template as<value_type>();
        penalbc = this->vm()["bccoeff"].template as<value_type>();
    }


    /**
     * create the mesh using mesh size \c meshSize
     */
    mesh_ptrtype createPerforedMesh( double meshSize );
    mesh_ptrtype createGlobalMesh( double meshSize );
    mesh_ptrtype createLocalMesh( double meshSize );


    /**
     * run the convergence test
     */
    void run();

private:

    template<typename IteratorRange, typename QuadratureT, typename ExprT,typename SpacePtrType>
    void integrateLocal( mesh_ptrtype localmesh,
                         mesh_ptrtype globalmesh,
                         SpacePtrType & Uh,
                         IteratorRange const& itrange, QuadratureT const& im, ExprT const& expr, vector_ptrtype& rhs );

    /**
     * solve system
     */
    void solve( sparse_matrix_ptrtype const& D, element_type& u, vector_ptrtype const& F, bool is_sym );
    void solveLocal( sparse_matrix_ptrtype const& D, element_type& u, vector_ptrtype const& F, bool is_sym );

    /**
     * export results to ensight format (enabled by  --export cmd line options)
     */
    void exportResults( element_type& u, uint time=0);
    void exportParticule( p0_element_type& U );
private:

    backend_ptrtype M_backend;
    backend_ptrtype M_backend_Loc;

    double meshSize,meshLocalSize;

    double mu;
    double penalbc;

    boost::shared_ptr<export_type> exporter;
}; // Stokes

template<int Dim, typename BasisU, typename BasisP, template<uint16_type,uint16_type,uint16_type> class Entity>
typename Stokes<Dim,BasisU,BasisP,Entity>::mesh_ptrtype
Stokes<Dim,BasisU,BasisP,Entity>::createPerforedMesh( double meshSize )
{
    mesh_ptrtype mesh( new mesh_type );

    GeoTool::Node x1(0,0);
    GeoTool::Node x2(1,1);
    GeoTool::Rectangle R( meshSize,"LEFT",x1,x2);
    R.setMarker(_type="line",_name="Lift",_marker3=true);
    R.setMarker(_type="line",_name="Paroi",_marker1=true,_marker2=true,_marker4=true);
    R.setMarker(_type="surface",_name="OmegaFluide",_marker1=true);

    GeoTool::Node x_centre(0.5,0.5);
    GeoTool::Node x_bord(0.5,0.7);
    GeoTool::Circle C( meshSize/6.,"disque",x_centre,x_bord);
    C.setMarker(_type="line",_name="Paroi",_markerAll=true);
    C.setMarker(_type="surface",_name="OmegaDisque",_marker1=true);


    //GmshTensorizedDomain<convex_type::nDim,convex_type::nOrder,convex_type::nRealDim,Entity> td;
    //td.setCharacteristicLength( meshSize );
    //std::string fname = td.generate( convex_type::name().c_str() );
    Gmsh gmsh;
    std::string fname = gmsh.generate( "domain", (R-C).geoStr()  );
    ImporterGmsh<mesh_type> import( fname );
    mesh->accept( import );
    return mesh;
} // Stokes::createMesh

template<int Dim, typename BasisU, typename BasisP, template<uint16_type,uint16_type,uint16_type> class Entity>
typename Stokes<Dim,BasisU,BasisP,Entity>::mesh_ptrtype
Stokes<Dim,BasisU,BasisP,Entity>::createGlobalMesh( double meshSize )
{
    mesh_ptrtype mesh( new mesh_type );

    GeoTool::Node x1(0,0);
    GeoTool::Node x2(1,1);
    GeoTool::Rectangle R( meshSize,"LEFT",x1,x2);
    R.setMarker(_type="line",_name="Lift",_marker3=true);
    R.setMarker(_type="line",_name="Paroi",_marker1=true,_marker2=true,_marker4=true);
    R.setMarker(_type="surface",_name="OmegaFluide",_marker1=true);

    Gmsh gmsh;
    std::string fname = gmsh.generate( "domain", R.geoStr()  );
    ImporterGmsh<mesh_type> import( fname );
    mesh->accept( import );
    return mesh;
} // Stokes::createMesh

template<int Dim, typename BasisU, typename BasisP, template<uint16_type,uint16_type,uint16_type> class Entity>
typename Stokes<Dim,BasisU,BasisP,Entity>::mesh_ptrtype
Stokes<Dim,BasisU,BasisP,Entity>::createLocalMesh( double meshSize )
{
    mesh_ptrtype mesh( new mesh_type );

    GeoTool::Node x_centre(0.5,0.5);
    GeoTool::Node x_bord(0.5,0.8);
    GeoTool::Circle C1( meshSize/6.,"disque1",x_centre,x_bord);
    C1.setMarker(_type="line",_name="outer",_markerAll=true);
    C1.setMarker(_type="surface",_name="OmegaDisque1",_marker1=true);

    //GeoTool::Node x_centre(0.5,0.5);
    GeoTool::Node x_bord2(0.5,0.7);
    GeoTool::Circle C2( meshSize/6.,"disque2",x_centre,x_bord2);
    C2.setMarker(_type="line",_name="Paroi",_markerAll=true);
    C2.setMarker(_type="surface",_name="OmegaDisque2",_marker1=true);


    //GmshTensorizedDomain<convex_type::nDim,convex_type::nOrder,convex_type::nRealDim,Entity> td;
    //td.setCharacteristicLength( meshSize );
    //std::string fname = td.generate( convex_type::name().c_str() );
    Gmsh gmsh;
    std::string fname = gmsh.generate( "domain", (C1-C2).geoStr()  );
    ImporterGmsh<mesh_type> import( fname );
    mesh->accept( import );
    return mesh;
} // Stokes::createMesh


template<int Dim, typename BasisU, typename BasisP, template<uint16_type,uint16_type,uint16_type> class Entity>
void
Stokes<Dim, BasisU, BasisP, Entity>::run()
{
    if ( this->vm().count( "help" ) )
        {
            std::cout << this->optionsDescription() << "\n";
            return;
        }

    using namespace Life::vf;

    this->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/" )
                            % this->about().appName()
                            % convex_type::name()
                            % BasisU::nOrder
                            % this->vm()["hsize"].template as<double>()
                            );

    //----------------------------------------------------------------//

    mesh_ptrtype meshGlobal = createGlobalMesh( meshSize );
    mesh_ptrtype meshLocal = createLocalMesh( meshLocalSize );
    mesh_ptrtype meshPerfored = createPerforedMesh( meshSize );

    //----------------------------------------------------------------//

    space_ptrtype XhGlob = space_type::New( meshGlobal );
    space_ptrtype XhLoc = space_type::New( meshLocal );
    space_ptrtype XhPerf = space_type::New( meshPerfored );

    p0_functionspace_ptrtype XhGlobP0 = p0_functionspace_type::New( meshGlobal );

    opinterp_u_ptrtype Ih_u = opinterp_u_ptrtype( new opinterp_u_type( XhGlob->template functionSpace<0>(),
                                                                       XhLoc->template functionSpace<0>(),
                                                                       M_backend ) );
    opinterp_p_ptrtype Ih_p = opinterp_p_ptrtype( new opinterp_p_type( XhGlob->template functionSpace<1>(),
                                                                       XhLoc->template functionSpace<1>(),
                                                                       M_backend ) );
    opinterp_l_ptrtype Ih_l = opinterp_l_ptrtype( new opinterp_l_type( XhGlob->template functionSpace<2>(),
                                                                       XhLoc->template functionSpace<2>(),
                                                                       M_backend ) );

    //----------------------------------------------------------------//

    element_type UGlob( XhGlob, "uGlob" );
    element_type VGlob( XhGlob, "vGlob" );
    element_0_type uGlob = UGlob.template element<0>();
    element_0_type vGlob = VGlob.template element<0>();
    element_1_type pGlob = UGlob.template element<1>();
    element_1_type qGlob = VGlob.template element<1>();
    element_2_type lambdaGlob = UGlob.template element<2>();
    element_2_type nuGlob = VGlob.template element<2>();

    //----------------------------------------------------------------//

    auto chip = chi( (Px()-cst(0.5))*(Px()-cst(0.5))+(Py()-cst(0.5))*(Py()-cst(0.5)) <= 0.1*0.1 );
    p0_element_type particule( XhGlobP0, "particule" );
    particule = vf::project( XhGlobP0, elements( meshGlobal ), chip );
    this->exportParticule( particule );

    // local problem ---------------------------------------------------------------------------
    element_type ULoc( XhLoc, "uLoc" );
    element_type VLoc( XhLoc, "vLoc" );
    element_0_type uLoc = ULoc.template element<0>();
    element_0_type vLoc = VLoc.template element<0>();
    element_1_type pLoc = ULoc.template element<1>();
    element_1_type qLoc = VLoc.template element<1>();
    element_2_type lambdaLoc = ULoc.template element<2>();
    element_2_type nuLoc = VLoc.template element<2>();

    auto deftLoc = gradt(uLoc);
    auto defLoc = grad(vLoc);
    auto defvLoc = gradv(uLoc);
    auto SigmaNtLoc = -idt(pLoc)*N()+mu*deftLoc*N();
    auto SigmaNLoc = -id(pLoc)*N()+mu*defLoc*N();
    auto SigmaNvLoc = -idv(pLoc)*N()+mu*defvLoc*N();

    sparse_matrix_ptrtype DLoc( M_backend->newMatrix( XhLoc, XhLoc ) );
    sparse_matrix_ptrtype DLocD( M_backend->newMatrix( XhLoc, XhLoc ) );

    // dummy matrix
    form2( XhLoc, XhLoc, DLoc, _init=true )=integrate( elements(meshLocal), mu*trace(deftLoc*trans(defLoc)) );
    form2( XhLoc, XhLoc, DLoc )+=integrate( elements(meshLocal), - div(vLoc)*idt(pLoc) + divt(uLoc)*id(qLoc) );
    form2( XhLoc, XhLoc, DLoc )+=integrate( elements(meshLocal), id(qLoc)*idt(lambdaLoc) + idt(pLoc)*id(nuLoc) );
    form2( XhLoc, XhLoc, DLoc )+=integrate( boundaryfaces(meshLocal),
                                   -trans(SigmaNtLoc)*id(vLoc)
                                   -trans(SigmaNLoc)*idt(uLoc)
                                   +penalbc*trans(idt(uLoc))*id(vLoc)/hFace() );

    DLoc->close();

    auto f = vec(cst(0),cst(0));
    auto vit=vec( Px()*-(1-Px()) , cst(0.));

    vector_ptrtype FLoc( M_backend->newVector( XhLoc ) );
    vector_ptrtype FLocD( M_backend->newVector( XhLoc ) );

    form1( XhLoc, FLoc, _init=true )  =  integrate( elements(meshLocal), trans(f)*id(vLoc) );

    FLoc->close();

   //----------------------------------------------------------------//

    auto deft = gradt(uGlob);
    auto def = grad(vGlob);
    auto defv = gradv(uGlob);
    auto SigmaNt = -idt(pGlob)*N()+mu*deft*N();
    auto SigmaN = -id(pGlob)*N()+mu*def*N();
    auto SigmaNv = -idv(pGlob)*N()+mu*defv*N();


    sparse_matrix_ptrtype DGlob( M_backend->newMatrix( XhGlob, XhGlob ) );
    sparse_matrix_ptrtype DGlobD( M_backend->newMatrix( XhGlob, XhGlob ) );

    form2( XhGlob, XhGlob, DGlob, _init=true )=integrate( elements(meshGlobal), mu*trace(deft*trans(def)) );
    form2( XhGlob, XhGlob, DGlob )+=integrate( elements(meshGlobal), - div(vGlob)*idt(pGlob) + divt(uGlob)*id(qGlob) );
    form2( XhGlob, XhGlob, DGlob )+=integrate( elements(meshGlobal), id(qGlob)*idt(lambdaGlob) + idt(pGlob)*id(nuGlob) );
    form2( XhGlob, XhGlob, DGlob)+=integrate( boundaryfaces(meshGlobal),
                                   -trans(SigmaNt)*id(vGlob)
                                   -trans(SigmaN)*idt(uGlob)
                                   +penalbc*trans(idt(uGlob))*id(vGlob)/hFace() );
    DGlob->close();

    vector_ptrtype FGlob( M_backend->newVector( XhGlob ) );
    vector_ptrtype FGlobD( M_backend->newVector( XhGlob ) );
    form1( XhGlob, FGlob, _init=true )  =
        integrate( elements(meshGlobal), trans(f)*id(vGlob) );
    form1( XhGlob, FGlob ) +=
        integrate( markedfaces(meshGlobal,meshGlobal->markerName("Lift")), trans(vit)*(-SigmaN+penalbc*id(vGlob)/hFace() ) );
    FGlob->close();

   //----------------------------------------------------------------//

#if 1
    // relaxation: note that using a LU solver should make it very
    // fast here at each iteration, only the right hand side is
    // actually updated
    int i = 1;
    double theta = 0.8;
    double eps = 1;
    double fixpointtol= 1.e-6;
    uint n_iter_relax = this->vm()["nbiter"].template as<int>();

    element_type ULocGamma( XhLoc, "varLoc" );
    element_0_type uLocGamma = ULocGamma.template element<0>();
    element_1_type pLocGamma = ULocGamma.template element<1>();
    element_2_type lambdaLocGamma = ULocGamma.template element<2>();

    std::cout<<"\nuSTART : "<<uLocGamma.start();
    std::cout<<"\npSTART : "<<pLocGamma.start();

    element_type ULocTemp( XhLoc, "varLoc" );
    element_0_type uLocTemp = ULocTemp.template element<0>();
    element_1_type pLocTemp = ULocTemp.template element<1>();
    element_2_type lambdaLocTemp = ULocTemp.template element<2>();
    ULocTemp.zero();
    ULocGamma.zero();


    element_type dnv( XhLoc, "u" );
    element_0_type dnv_u = dnv.template element<0>();
    element_type dnv1( XhLoc, "u" );
    element_0_type dnv1_u = dnv1.template element<0>();
    element_type dnvm1( XhLoc, "u" );
    element_0_type dnvm1_u = dnvm1.template element<0>();

    dnv.zero();
    dnv1.zero();
    dnvm1.zero();

    do
    {
        boost::timer ti;
        std::cout << "============================================================\n";
        std::cout << "[FBM] iteration = " << i << "\n";

        //relaxation : ULoc = theta*ULoc0 + (1-theta)*u;
        ULocTemp.zero();
        ULocTemp.add(1,ULocGamma);
        ULocGamma.zero();
        ULocGamma.add( theta, ULocTemp );

        Ih_u->apply( uGlob, uLocTemp );
        Ih_p->apply( pGlob, pLocTemp );
        Ih_l->apply( lambdaGlob, lambdaLocTemp );

        ULocGamma.add( 1-theta, ULocTemp );
        std::cout << "relaxation done in " << ti.elapsed() << "s\n";ti.restart();

        ///////////////
        DLocD->close();
        DLocD->zero();
        DLocD->addMatrix(1.0, DLoc );

        FLocD->close();
        FLocD->zero();
        FLocD->add(1.0, FLoc );

        if (true )// this->vm().count( "weak" ) )
        {
            form1( XhLoc, FLocD ) +=
                integrate( markedfaces( meshLocal, meshLocal->markerName( "outer" ) ), // _Q<max20<4*Order>::order>(),
                           trans(idv(uLocGamma) )*(-SigmaNLoc+penalbc*id(vLoc)/hFace() ) );
        }

        DLocD->close();
        FLocD->close();
#if 0
        if (true)// !this->vm().count( "weak" ) )
        {
            form2( XhLoc, XhLoc, DLocD ) += on( markedfaces( meshLocal, meshLocal->markerName( "Lift" ) ), uLoc, FLoc, idv(uLocGamma) );
            form2( XhLoc, XhLoc, DLocD ) += on( markedfaces( meshLocal, meshLocal->markerName( "Paroi" ) ), uLoc, FLoc, constant(0.) );
        }
#endif

        std::cout << "local rhs done in " << ti.elapsed() << "s\n";ti.restart();

        this->solveLocal( DLocD, ULoc, FLocD, false );
        std::cout << "sove local system done in " << ti.elapsed() << "s\n";ti.restart();


        if ( i == 1 )
            dnv1 = ULoc;
        dnvm1=dnv;
        dnv = ULoc;



        FGlobD->close();
        FGlobD->zero();
        FGlobD->add( 1.0, FGlob );

        //auto oio = space_p_type::New(meshGlobal);
        auto oio = space_u_type::New(meshGlobal);

        //typedef typename _Q<max20<4*Order>::order>::template apply<mesh_type::nDim,typename mesh_type::value_type, Simplex>::type quad_type2;
        typedef typename _Q<6>::template apply<mesh_type::nDim,typename mesh_type::value_type, Simplex>::type quad_type2;
        quad_type2 __quadra;
        this->integrateLocal//<space_p_type>
            ( meshLocal,
              meshGlobal,
              oio,//pGlob.functionSpace(),//XhGlob->template functionSpace<1>(),
              markedfaces( meshLocal, meshLocal->markerName( /*"inner"*/"Paroi" )),
              __quadra,
              //grad_g*N(), FGlobFull );
              //grad_gGamma*N(), FGlobFull );
              //gradv(pGlob)*N(), FGlobD );
              //gradv(uGlob)*N(), FGlobD);
              SigmaNvLoc, FGlobD );


        //std::cout << "\ninteg sigma : " << integrate(markedfaces( meshLocal, meshLocal->markerName( "Paroi" )),
        //                                             -trans(idv(pLoc)*N() + mu*gradv(uLoc)*N())*(idv(pLoc)*N() + mu*gradv(uLoc)*N())   ).evaluate()(0,0);
        /*
        std::cout << "\ninteg sigma : " << integrate(markedfaces( meshLocal, meshLocal->markerName( "Paroi" )),
                                                     -trans(idv(pLoc)*N())*idv(pLoc)*N()
                                                     + mu*trans(gradv(uGlob)*N())*(gradv(uGlob)*N())  ).evaluate()(0,0);
        */
        //std::cout << "\ninteg sigma : " << integrate(elements( meshLocal),trans(SigmaNvLoc)*SigmaNvLoc ).evaluate()(0,0);
        //std::cout << "\ninteg sigma : " << integrate(markedfaces( meshLocal, meshLocal->markerName( "Paroi" )),trans(SigmaNvLoc)*SigmaNvLoc ).evaluate()(0,0);
        std::cout <<"\nNorme FGlobD : "<< FGlobD->l2Norm() <<"\n";

        DGlobD->close();
        DGlobD->zero();
        DGlobD->addMatrix(1.0, DGlob );

        this->solve( DGlobD, UGlob, FGlobD, false );

        std::cout << "global problem solve done in " << ti.elapsed() << "s\n";ti.restart();

        //std::cout << "export done in " << ti.elapsed() << "s\n";ti.restart();
        // next step
        //eps = std::sqrt(integrate( markedfaces( meshLocal, meshLocal->markerName( "Paroi" )) ,// _Q<6>(),
            //                         trans(gradv(dnv_u)*N()-gradv(dnvm1_u)*N() )*(gradv(dnv_u)*N()-gradv(dnvm1_u)*N()) ).evaluate()(0,0)/* /int_gradgN*/);

        if ( i > 1 )
            eps = std::abs(integrate( markedfaces( meshLocal, meshLocal->markerName( "Paroi" )) ,// _Q<6>(),
                                      idv(dnv_u)*N()-idv(dnvm1_u)*N()
                                      + gradv(dnv_u)*N()-gradv(dnvm1_u)*N() ).evaluate()(0,0)/* /int_gradgN*/);

        std::cout << "fix point error = " << eps << "\n";
        this->exportResults( UGlob,i-1);
        this->exportResults( ULoc,i-1);
        ++i;

    }
    while( (eps > fixpointtol ) && (i < n_iter_relax ) );
#endif

    //this->solve( DGlob, UGlob, FGlob, false );

    //this->exportResults( UGlob);
    //this->exportResults( ULoc);

} // Stokes::run

template<int Dim, typename BasisU, typename BasisP, template<uint16_type,uint16_type,uint16_type> class Entity>
void
Stokes<Dim, BasisU, BasisP, Entity>::solve( sparse_matrix_ptrtype const& D,
                                   element_type& u,
                                   vector_ptrtype const& F,
                                   bool is_sym )
{
    vector_ptrtype U( M_backend->newVector( u.functionSpace() ) );
    M_backend->solve( D, D, U, F, false );
    u = *U;
} // Stokes::solve


template<int Dim, typename BasisU, typename BasisP, template<uint16_type,uint16_type,uint16_type> class Entity>
void
Stokes<Dim, BasisU, BasisP, Entity>::solveLocal( sparse_matrix_ptrtype const& D,
                                                 element_type& u,
                                                 vector_ptrtype const& F,
                                                 bool is_sym )
{
    vector_ptrtype U( M_backend_Loc->newVector( u.functionSpace() ) );
    M_backend_Loc->solve( D, D, U, F );
    u = *U;
} // Fbm::solve


template<int Dim, typename BasisU, typename BasisP, template<uint16_type,uint16_type,uint16_type> class Entity>
void
Stokes<Dim, BasisU, BasisP, Entity>::exportResults( element_type& U,uint time )
{

    export_ptrtype __exporter( export_type::New( "ensight",//"gmsh",
                                                 this->about().appName()+U.name() ) );
    //__exporter->step( 0 )->setMesh( U.template element<0>().functionSpace()->mesh() );
    __exporter->step( time )->setMesh( U.functionSpace()->mesh() );
    __exporter->step( time )->add( "u", U.template element<0>() );
    __exporter->step( time )->add( "p", U.template element<1>() );
    __exporter->save();

} // Stokes::export
template<int Dim, typename BasisU, typename BasisP, template<uint16_type,uint16_type,uint16_type> class Entity>
void
Stokes<Dim, BasisU, BasisP, Entity>::exportParticule( p0_element_type& U )
{
    export_ptrtype __exporter( export_type::New( "ensight",//"gmsh",
                                                 this->about().appName()+"particule" ) );

    __exporter->step( 0 )->setMesh( U.functionSpace()->mesh() );
    __exporter->step( 0 )->add( "particule", U );
    __exporter->save();

} // Stokes::export

template<int Dim, typename BasisU, typename BasisP, template<uint16_type,uint16_type,uint16_type> class Entity>
template<typename IteratorRange, typename QuadratureT, typename ExprT,typename SpacePtrType>
void
Stokes<Dim, BasisU, BasisP, Entity>::integrateLocal( mesh_ptrtype localmesh,
                                                     mesh_ptrtype globalmesh,
                                                     SpacePtrType & Uh,
                                                     IteratorRange const& itrange, QuadratureT const& im, ExprT const& expr, vector_ptrtype& rhs )
{
    //
    // some typedefs
    //
    typedef typename mesh_type::element_type geoelement_type;
    typedef typename mesh_type::element_iterator mesh_element_iterator;
    typedef typename mesh_type::gm_type gm_type;

    typedef mesh_type local_mesh_type;
    typedef typename local_mesh_type::element_type local_geoelement_type;
    typedef typename local_mesh_type::element_iterator local_mesh_element_iterator;
    typedef typename local_mesh_type::gm_type local_gm_type;

    // local
    typedef typename local_gm_type::template Context<ExprT::context|vm::JACOBIAN|vm::KB|vm::NORMAL|vm::POINT, local_geoelement_type> local_gmc_type;
    typedef boost::shared_ptr<local_gmc_type> local_gmc_ptrtype;

    // global
    typedef typename gm_type::template Context<vm::POINT, geoelement_type> global_gmc_type;
    typedef boost::shared_ptr<global_gmc_type> global_gmc_ptrtype;

    //
    // Precompute some data in the reference element for
    // geometric mapping and reference finite element
    //
    typedef typename gm_type::precompute_type pc_type;
    typedef typename gm_type::precompute_ptrtype pc_ptrtype;

    typedef typename local_gm_type::precompute_type local_pc_type;
    typedef typename local_gm_type::precompute_ptrtype local_pc_ptrtype;



    typedef typename SpacePtrType::element_type space_ctx_type;


    //typedef typename QuadratureT::template apply<mesh_type::nDim,typename mesh_type::value_type, Simplex>::type quad_type;

    //    QuadMapped<quad_type> qm;
    //typedef typename QuadMapped<quad_type>::permutation_type permutation_type;
    //typename QuadMapped<quad_type>::permutation_points_type ppts( qm( im ) );

    QuadMapped<QuadratureT> qm;
    typedef typename QuadMapped<QuadratureT>::permutation_type permutation_type;
    typename QuadMapped<QuadratureT>::permutation_points_type ppts( qm( im ) );

    std::vector<std::map<permutation_type, local_pc_ptrtype> > __geopc( im.nFaces() );


    //    typedef typename quad_type::face_quadrature_type face_im_type;
    typedef typename QuadratureT::face_quadrature_type face_im_type;
    //typedef typename QuadratureT::template apply<mesh_type::nDim,typename mesh_type::value_type, Simplex>::type::face_quadrature_type face_im_type;

    //__typeof__(im(__face_id_in_elt_0 ) ) im_face ( im(__face_id_in_elt_0 ) );
    std::vector<face_im_type> face_ims( im.nFaces() );

    for ( uint16_type __f = 0; __f < im.nFaces(); ++__f )
    {
        face_ims[__f] = im.face( __f );

        for( permutation_type __p( permutation_type::IDENTITY );
             __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
        {
            __geopc[__f][__p] = local_pc_ptrtype(  new local_pc_type( localmesh->gm(), ppts[__f].find(__p)->second ) );
        }
    }

    Log() << "[FBM] integrate Local\n";

    typedef typename boost::tuples::template element<1, IteratorRange>::type iterator_type;
    iterator_type local_it = itrange.template get<1>();
    iterator_type local_en = itrange.template get<2>();

    Log() << "[FBM::integrateLocal] n elts: " << std::distance( local_it, local_en ) << "\n";
    uint16_type __face_id_in_elt_0 = local_it->pos_first();

    // check that we have elements to iterate over
    if ( local_it == local_en )
        return;

    // get the geometric mapping associated with element 0
    local_gmc_ptrtype local_c0( new local_gmc_type( localmesh->gm(),
                                                    local_it->element( 0 ),
                                                    __geopc,//__geopc[__face_id_in_elt_0][local_it->element(0).permutation( __face_id_in_elt_0 )],
                                                    __face_id_in_elt_0 ) );

    uint16_type nptsface = im.nPointsOnFace();

    // locate the points in the global mesh
    typename mesh_type::Inverse meshinv( globalmesh );

    size_type nelts = std::distance(local_it,local_en);
    //std::vector<double> expr_value( nelts*nptsface );
    typedef boost::multi_array<double,2> array_type;
    array_type expr_value( boost::extents[space_ctx_type::nComponents1][nelts*nptsface] );

    for( size_type ide = 0; local_it != local_en; ++ local_it, ++ide )
    {
        //Log() << "[FBM] local element " << local_it->id() << "\n";
        // generate the points in the real space
        uint16_type __face_id_in_elt_0 = local_it->pos_first();
        local_c0->update( local_it->element(0),
                          //__geopc[__face_id_in_elt_0][local_it->element(0).permutation(__face_id_in_elt_0)],
                          __face_id_in_elt_0 );

        typedef fusion::map<fusion::pair<vf::detail::gmc<0>, local_gmc_ptrtype> > map_gmc_type;

        typedef typename ExprT::template tensor<map_gmc_type> eval_expr_type;
        typedef boost::shared_ptr<eval_expr_type> eval_expr_ptrtype;
        map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( local_c0 ) );
        eval_expr_ptrtype eval_expr( new eval_expr_type( expr, mapgmc ) );


        eval_expr->init( im );
		eval_expr->update( mapgmc, __face_id_in_elt_0 );
        for( int q = 0; q < nptsface; ++ q )
        {
            size_type idq = nptsface*ide+q;
            // evaluate expression and store it
            for (uint c=0;c<space_ctx_type::nComponents1;++c)
                expr_value[c][idq] = face_ims[__face_id_in_elt_0].weight(q)*eval_expr->evalq( 0, 0, q )*local_c0->J(q)*local_c0->normalNorm(q);
            // insert quadrature point in inverse mesh data structure
            meshinv.addPointWithId( local_c0->xReal( q ), idq, 0 );
        }
    }
    meshinv.distribute();



    std::vector<bool> dof_done( nelts*nptsface );
    std::fill( dof_done.begin(), dof_done.end(), false );
    std::vector<boost::tuple<size_type,uint16_type > > itab;
    matrix_node<value_type>::type pts( mesh_type::nDim, 1 );


    mesh_element_iterator globalit = globalmesh->beginElementWithProcessId( Application::processId() );
    mesh_element_iterator globalen = globalmesh->endElementWithProcessId( Application::processId() );

    

    // basis
    //typedef typename space_type::basis_type basis_type;
    typedef typename space_ctx_type::basis_type basis_type;
    //typedef typename space_type::template sub_functionspace<0>::type::element_type::basis_type basis_type;
    basis_type const* global_basis = Uh->basis().get();
    pc_ptrtype global_geopc( new pc_type( Uh->gm(), global_basis->dual().points() ) );
    global_gmc_ptrtype global_gmc( new global_gmc_type( Uh->gm(), *globalit, global_geopc ) );
    typedef typename basis_type::template Context<global_gmc_type::context,basis_type,gm_type,geoelement_type> basiscontext_type;
    typedef boost::shared_ptr<basiscontext_type> basiscontext_ptrtype;


    //!!!!! _M_Xh->nDofStart( _M_index ), 0 ) );

    for( ; globalit != globalen; ++ globalit )
    {
        global_gmc->update( *globalit/*, global_geopc*/ );
        meshinv.pointsInConvex( globalit->id(), itab );
        if (itab.size() == 0)
            continue;

        for (size_type i = 0; i < itab.size(); ++i)
        {
            // get dof id in target dof table
            size_type q = boost::get<0>( itab[i] );
            if ( !dof_done[q] )
            {
                dof_done[q]=true;
                ublas::column( pts, 0 ) = meshinv.referenceCoords()[q];

                global_geopc->update( pts );
                global_gmc_ptrtype global_gmccc( new global_gmc_type( Uh->gm(), *globalit, global_geopc ));
                //global_gmc->update( *globalit/*, global_geopc*/ );//Y a un souci?

                //typedef typename element_type::pc_type elementpc_type;
                //typedef typename element_type::pc_ptrtype elementpc_ptrtype;
                typedef typename space_ctx_type::element_type::pc_type elementpc_type;
                typedef typename space_ctx_type::element_type::pc_ptrtype elementpc_ptrtype;

                elementpc_ptrtype pc( new elementpc_type( Uh->fe(), pts ) );
                basiscontext_ptrtype global_bc( new basiscontext_type( Uh->fe(), global_gmccc, pc ) );

                std::cout<<"bas"<<basis_type::nComponents;
                for (uint16_type loc_ind=0; loc_ind < basis_type::nLocalDof ; loc_ind++)
                {
                    for ( uint16_type comp = 0;comp < basis_type::nComponents;++comp )
                    {
                        size_type globaldof = boost::get<0>(Uh->dof()->localToGlobal(globalit->id(), loc_ind, comp ));

                        rhs->add( globaldof, expr_value[comp][q]*global_bc->id( loc_ind, 0, comp, 0 ) );
                    }
                }

            }
        }
    } // global_it

    // make sure that all quad points have been visited
    for( size_type q = 0; q < dof_done.size(); ++q )
    {
        if ( dof_done[q] != true )
        {
            Log() << "quad point " << q << " has not been found in global mesh\n";
        }
    }
    return;
    //#endif
}




} // Life
} // Cavity
int
main( int argc, char** argv )
{

    using namespace Life;
    /* assertions handling */
    Life::Assert::setLog( "stokes.assert");

    const int nDim = 2;
    typedef FBM_Cavity::Stokes<nDim, Lagrange<2, Vectorial>,Lagrange<1, Scalar>, Simplex> stokes_type;
    //typedef Life::Stokes<nDim, CrouzeixRaviart<1, Vectorial>,Lagrange<0, Scalar>, Simplex> stokes_type;


    /* define and run application */
    stokes_type stokes( argc, argv, FBM_Cavity::makeAbout(), FBM_Cavity::makeOptions() );
    stokes.run();
}





